#!/usr/bin/perl
# 2014/01/07 just BWA
# 2012/09/26 just convert BWA output file into bowtie format


use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;


#---------------------------------------
# length of each reads
#---------------------------------------
my $LENGTH = 50;


my %param = (
	2	=> 1,
	4	=> 1,
);

if( !(exists $param{scalar(@ARGV)}) or $ARGV[0] eq '--help'){
	my $help =<<"EOF";
Usage : $0 -i [sanger fastq file] -j [jump position]
EOF
	die "$help";
}

my %opt;
getopts("i:j:", \%opt);
my $FILE_fastq = $opt{i};
(my $FILE_sai = $FILE_fastq) =~ s/\.fastq/.sai/;
(my $FILE_sam = $FILE_fastq) =~ s/\.fastq/.sam/;
(my $FILE_bowtie = $FILE_fastq) =~ s/\.fastq/_bowtie.map/;
(my $FILE_alnLog = $FILE_fastq) =~ s/\.fastq/_aln.log/;
(my $FILE_samLog = $FILE_fastq) =~ s/\.fastq/_sam.log/;
my $label = $opt{j};

#my $FILE_reference = "/usr/local/biology/bwa-0.6.2/pombe2009";
my $FILE_reference = "/home/hideki/bio/bwa/pombe.ASM294v2.19";

### jump to some position
if(exists $opt{j}){
	goto $label;
}


#---------------------------------------
# aln
#---------------------------------------
system("bwa aln -t 4 $FILE_reference $FILE_fastq > $FILE_sai 2> $FILE_alnLog") == 0 or die "aln command failed: $!";


#---------------------------------------
# samse
#---------------------------------------
system("bwa samse -n 100 $FILE_reference $FILE_sai $FILE_fastq > $FILE_sam 2> $FILE_samLog") == 0 or die "samse command failed: $!";

CONVERT:
#---------------------------------------
# SAM format
#---------------------------------------
#@SQ	SN:chrI	LN:5579133
#@SQ	SN:chrII	LN:4539804
#@SQ	SN:chrIII	LN:2452883
#@SQ	SN:chrM	LN:19431
#@PG	ID:bwa	PN:bwa	VN:0.5.9-r16
#ILLUMINA-E6BD77_0034_FC:3:1:1749:1125#0	4	*	0	0	*	*	0	0	GACTGNTGAATTCAAACAATGGCGAGAAAAGATATG	BB==7#?9DD?BEDEGGGEEGGFGGG=GDDA@DEGG
#ILLUMINA-E6BD77_0034_FC:3:1:2292:1128#0	16	chrI	4235288	25	36M	*	0	0	ATAGCATCGTATAGTATCGTATCGTATCGCNATAGA	@DBDDD2@DDDD@DDBDBDDB8BB=B?BBB#?A@@A	XT:A:U	NM:i:2	X0:i:1	X1:i:0	XM:i:2	XO:i:0	XG:i:0	MD:Z:12C17C5
#ILLUMINA-E6BD77_0034_FC:3:1:2416:1129#0	0	chrI	3074849	37	36M	*	0	0	TGTCTNGTGTTTAAGCCCTTCAATATGCGTGTAAAG	CBCCC#FEFFFFFFFHHHHHHFHHHHHHHHDBGGEH	XT:A:U	NM:i:1	X0:i:1	X1:i:0	XM:i:1	XO:i:0	XG:i:0	MD:Z:5T30
#ILLUMINA-E6BD77_0034_FC:3:1:2683:1124#0	0	chrI	2908596	0	36M	*	0	0	AAATCNTTCAAGCAAATCTCCAAGTGGAGCTCACCA	EEEDE#GGGGGGGGGIIIIIIIIGIIIHGIIIIIIH	XT:A:R	NM:i:1	X0:i:2	X1:i:0	XM:i:1	XO:i:0	XG:i:0	MD:Z:5C30	XA:Z:chrIII,-538924,36M,1;
#ILLUMINA-E6BD77_0034_FC:3:1:2752:1128#0	16	chrII	3984601	37	36M	*	0	0	TCTTCCGTTTCGTCTGTTTTTCCTGTTAACNATTCC	GDIIHGIII@IIIHIIIIIIIGGGGGFFEF#DEDEE	XT:A:U	NM:i:1	X0:i:1	X1:i:0	XM:i:1	XO:i:0	XG:i:0	MD:Z:30T5
#ILLUMINA-E6BD77_0034_FC:3:1:3052:1125#0	0	chrIII	1229144	37	36M	*	0	0	GATTANTACCTTTTTAATTTAATTTTGTCCTACCTT	C?CCA#FFEFGGGGGGIIIIIIIIIIGIIIIIIIII	XT:A:U	NM:i:1	X0:i:1	X1:i:0	XM:i:1	XO:i:0	XG:i:0	MD:Z:5A30
#ILLUMINA-E6BD77_0034_FC:3:1:3274:1126#0	4	*	0	0	*	*	0	0	AAAGANGATGAATGGATTAAAAGGTGTTGGAGTAGA	EEDEC#FEFFEEEFFGIIIIGGGGFGGGGGIIIGGG


#---------------------------------------
# fake Bowtie format 最後のカラムはリピートの数
#---------------------------------------
#ILLUMINA-E6BD77_0012:3:1:1377:20538#0/1 -       chr3    1901683 1
#ILLUMINA-E6BD77_0012:3:1:1376:17697#0/1 -       chr2    167853  1
#ILLUMINA-E6BD77_0012:3:1:1377:18765#0/1 +       chr1    727014  1
#ILLUMINA-E6BD77_0012:3:1:1377:14335#0/1 +       chr2    492656  1
#ILLUMINA-E6BD77_0012:3:1:1377:13143#0/1 -       chr2    4417538 13

my $TOTAL_read = 0;
my $Aligned_read = 0;

my $fh_sam = IO::File->new($FILE_sam) or die "cannot open $FILE_sam: $!";
my $fh_out = IO::File->new($FILE_bowtie, 'w') or die "cannot write $FILE_bowtie: $!";
while($_ = $fh_sam->getline()){
	s/\r?\n//;
	if(m/^\@/){
		next;
	}
	my ($id, $flag, $chr, $position, $mapQ, $CIAGR, $mate_a, $mate_b, $mate_c, $seq, $quality, @option) = split /\t/;
	my ($xt, $x0, $xa) = &CheckOption(@option);

	$TOTAL_read++;
	if($flag == 0){
		my $direction = '+';
		$fh_out->printf("%s\t%s\t%s\t%d\t%d\n", $id, $direction, $chr, $position, $x0);
		$Aligned_read++;
	}elsif($flag == 16){
		my $direction = '-';
		my $bowtie_position = $position + $LENGTH - 1;
		$fh_out->printf("%s\t%s\t%s\t%d\t%d\n", $id, $direction, $chr, $bowtie_position, $x0);
		$Aligned_read++;
	}else{
		# flag == 4 (unmapped)
		next;
	}
}
$fh_sam->close();
$fh_out->close();

printf "%d / %d = %.2f\n", $Aligned_read, $TOTAL_read, $Aligned_read / $TOTAL_read;


sub CheckOption{
	my @data = @_;
	my $xt = 'X';
	my $x0 = 0;
	my $xa = 'NA';
	foreach my $d(@data){
		if($d =~ m/^XT:A:(\w)/){
			$xt = $1;
		}elsif($d =~ m/^X0:i:(\d{1,3})/){
			$x0 = $1;
		}elsif($d =~ m/^XA:Z:(.+)/){
			$xa = $1;
		}
	}
	return ($xt, $x0, $xa);
}


